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Frequency assortativity can induce chaos in oscillator networks 
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We investigate the effect of preferentially connecting oscillators with similar frequency to each other in net¬ 
works of coupled phase oscillators (i.e., frequency assortativity). Using the network Kuramoto model as an 
example, we find that frequency assortativity can induce chaos in the macroscopic dynamics. By applying 
a mean-field approximation in combination with the dimension reduction method of Ott and Antonsen, we 
show that the dynamics can be described by a low dimensional system of equations. We use the reduced sys¬ 
tem to characterize the macroscopic chaos using Lyapunov exponents, bifurcation diagrams, and time-delay 
embeddings. Finally, we show that the emergence of chaos stems from the formation of multiple groups of 
synchronized oscillators, i.e., meta-oscillators. 

PACS numbers: 05.45.Xt, 89.75.Hc 


The synchronization of network-coupled dynamical sys¬ 
tems |[Il,|2l plays a key role in many natural phenomena |l3l01 
and engineering applications nil]. An important example is 
networks of coupled oscillators. Kuramoto showed Jj] that 
under suitable conditions, the analysis of an ensemble of N 
oscillators can be reduced to the dynamics of phase angles 
for the oscillators, where oscillator i has phase angle 6i for 
z = 1,..., W. When the oscillators are coupled by a network, 
the corresponding model is given by 

N 

= Wi + a: ^ Ay sin (6»y - 0i), (1) 

where Wi is the natural frequency of oscillator i, K > 0 is the 
global coupling strength, and [Ay] is the network adjacency 
matrix that encodes the network structure (Ay = 1 if there is 
a network link from node j to node i and Ay = 0 otherwise). 

The dynamics of Eq. O and its many extensions have 
been the subject of a great deal of research (e.g., Refs. IH- 
d). Recently an advance in the analysis of such systems 
was obtained IUaIl^ which posits an ansatz for the long time 
asymptotic form of the solution of such systems and results in 
a dimensionality reduction whereby the W-dimensional dy¬ 
namics of Eq. O can be reduced to a much smaller system. 
This ansatz was first used on all-to-all coupled phase oscil¬ 
lator systems 0 (where each entry of the adjacency matrix 
is Ay = 1), and adapted to obtain analytical results reveal¬ 
ing the effects of various extensions of the original Kuramoto 
model, including chimera states, periodic forcing, bimodal 
frequency distributions, time-delays, clustering, and commu¬ 
nities 111514^ . Recently, the ansatz was extended via a mean- 
field technique to allow for the treatment of nontrivial network 
topologies d^, importantly shedding light on the effects of 
correlations between the degrees of network-connected node 
pairs, i.e., degree assortativity ll^ . 
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The formalism of Ref. can in principle be extended to 
account for assortativity based on arbitrary nodal properties, 
i.e., for probabilistic network generative models in which the 
probability that two nodes are connected is a function of pre¬ 
assigned nodal properties HQ. In particular, referring to 
Eq. ([T]i we note that nodes are characterized not only by their 
in- and out-degrees (fc™ = Ay, = Yj bnt 

also by their natural frequencies uji. It would seem that fre¬ 
quency assortativity would be crucial for the dynamics of the 
network Kuramoto problem since cooperative interactions be¬ 
tween pairs of connected nodes with like (unlike) frequencies 
would be stronger (weaker). However, so far there is no ana¬ 
lytical means of investigating the impact of this basic consid¬ 
eration on network dynamics. It is the purpose of this Rapid 
Communication to provide and illustrate such an analytical 
technique for investigating this effect. Our results show that 
frequency assortativity can play a profound role in determin¬ 
ing dynamical behavior. In particular, we show that frequency 
assortativity can induce chaos in the macroscopic system dy¬ 
namics. While chaos has previously been found in the macro¬ 
scopic dynamics of phase oscillator models we find it 

remarkable that chaos and complex dynamics can arise in the 
simple, basic model given by Eq. (|TJ merely from frequency 
assortativity. In the remainder of this Rapid Communication 
we describe a simple model for generating networks with fre¬ 
quency assortativity, investigate the emergence of chaotic dy¬ 
namics in such networks using numerical simulations, present 
a dimensionality reduction method for such networks, and fi¬ 
nally close with a brief discussion of our results. 

Frequency assortativity network model. We begin by 
briefly describing a model for generating oscillator net¬ 
works with particular frequency-frequency correlation be¬ 
tween neighbors, i.e., frequency assortativity. In other words, 
this model will allow for the construction of networks where 
neighboring oscillator tend to have similar or dissimilar nat¬ 
ural frequencies. Because we wish to focus on the effect of 
frequency assortativity in the simplest and cleanest setting, 
we henceforth consider the case of an undirected network in 
which all nodes have the same degree. Note that, by this 
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FIG. 1. (Color online) Synchronization in non-assortative and as- 
sortative networks. Synchronization profiles R (solid blue) and B 
(dashed red) vs K for examples of (a) non-assortative and (b) assor- 
tative networks of size N = 1000 with constant degree k = 50. For 
the reduced description on which our determination of B is based we 
use N = 20 and tUmax, —tUmin = 3.126. 


choice, issues of different degree distibutions, node degree- 
frequency correlations, and degree assortativity are, by defi¬ 
nition, absent, thus providing an unambiguous testing ground 
for investigating frequency assortativity effects with no other 
complications. (We note that, although our subsequenct the¬ 
ory is for this special case, it is easily generalized to ac¬ 
count for the other effects mentioned above.) Our model is 
based on the configuration model such that to each node 
i = 1,..., N we assign the same degree ki = k. Additionally, 
we assign to each oscillator a target frequency, wo.i which will 
be used to build network connections as follows. Choosing a 
node i that still requires at least one additional link, another 
node i which still requires at least one link is chosen accord¬ 
ing to a probability pij. Each pij depends on the target fre¬ 
quencies and Wqj- In the networks used here, we use 
Pij oc 0.5 + cld?/{(P + 1^0,i — wo,y|^) — 0.5] with d = 0.8 
and 7 = 5. In essence, the parameter c tunes the degree of fre¬ 
quency assortativity: c > 0 (c < 0) allows oscillators to more 
likely make connections to other oscillators with similar (dis¬ 
similar) target frequencies, resulting in assortative (disassorta- 
tive) networks. Links are made until all nodes have degree k. 
Finally, actual natural frequencies are assigned to each oscil¬ 
lator i according to a distribution . (w) that depends on the 
target frequency uoy. Here we consider the case of Lorentzian 
distributions 


where Ri{t) = describes the local order pa¬ 

rameter for oscillator i, we plot the evolution of R K in 
Fig-Efor the non-assortative and assortative networks in pan¬ 
els (a) and (b), respectively, using a solid blue curve. While 
the non-assortative networks displays typical behavior, transi¬ 
tioning from incoherence (R Ri 0) to coherence (R > 0) at a 
finite coupling strength (K Ki 0.05), the assortative network 
displays much more interesting behavior. In particular, in a 
range of intermediate coupling strengths (0.04 ^ K < 0.12) 
the order parameter undergoes large, irregular oscillations. 
We will now present a dimensionality reduction method which 
we will use to show that the dynamics in this interesting 
regime are in fact chaotic. 

Dimensionality reduction. The analytical technique we 
now summarize represents an extension of that described in 
Ref. Here we assume that the general structure de¬ 

scribed above, i.e., a network described by a single degree 
and a collection of target frequencies, the latter specifying 
the distributions . (w) from which the natural frequencies 
are drawn. The network is therefore characterized by the tar¬ 
get frequency distribution , which is normalized such that 
Swo ^<^0 — The frequency assortativity of the network is 
captured by the function the probability that a link 

exists from an oscillator with target frequency Wq to one with 
ujQ. We note that the assortativity function is con¬ 

strained to satisfy 

(4) 


We proceed by considering the limit of large networks, i.e., 
N ^ oo, such that the state of the network can be described 
by the family of distribution functions f^jg {6, uj, t), where 
fgjg{9,uj,t)d9duj/2Tr is the fraction of oscillators with target 
frequency ojq with phase in [9,9 -\- d0] and natural frequency 
in [w, a; + dw] at time t. We emphasize that each natural fre¬ 
quency depends on the target frequency, and since uj does not 
change in time we have 

^27r 

=p<^o(w). (5) 


9ujo{^) 


_ 

TT (w - UoY + A2^ ’ 


( 2 ) 


centered at wq with spread . 

We next demonstrate the effect of frequency assortativity by 
presenting results from numerical simulations. Considering a 
network of size N = 1000 with constant degree k = 50, we 
generate a non-assortative network and an assortative network 
using c = 0 and c = 1, respectively, and set A^^g = 0.05. 
We next solve Eq. O for each network, increasing K from 
zero by an increment of 10“® at each timestep At = 0.002. 
Defining the order parameter 


R{t) 


1 


N 


(3) 


The interaction term in Eq. O for an oscillator j can 
be expressed in terms of the local order parameters as 
ATm(e“*®^ Rj). The mean-field version of the local order pa¬ 
rameter is Ri{t) Rgjg ^ (f) and is given by 

Rojoii) = '^Puj'gau,’g^u,o Jj ( 6 ) 

Ul'g 

Finally, by the conservation of the number of oscillators, each 
distribution must satisfy the continuity equation 

0 = dtfujoi9,uj,t) + dg[{uj + Rgjg{t)]) 

(7) 

Together, Eqs. @ and (|7]) give a mean-field description for the 
macroscopic dynamics of Eq. O- (We note that in other con- 



























3 


texts the assortativity can be formulated in terms of degrees 
by replacing with IE^Ij or, still more gener¬ 

ally, ak',u; Q ■) 

We now follow Refs. EH where the authors showed 
that in the long-time limit each distribution function ap¬ 
proaches the form 


f^g{0,UJ,t) 




n—1 


( 8 ) 


where c.c. denotes the complex conjugate of the preceding 
term. [Note that, since ® is the time asymptotic form of the 
distribution, our use of ® should yield a good approximation 
of all the attractor dynamics, but not necessarily the transient 
dynamics that describes the approach to an attractor.] Substi¬ 
tuting Eq. (O in Eq. (|7]), we find that each satisfies 

dtb^^ {uj,t)= iujb^„ (w, f) -f y (f) - bl^ (w, (f)] . 

(9) 

Next, we substitute Eq. (|8ll into Eq. (|6]l to obtain 

Rcooit) = J 9uj'„{uj)buj'^{u}' ,t)duj'. ( 10 ) 

Assuming that the frequency distribution are Lorentzian as 
in Eq. dUl, Eg. (ITOl i can be simplified using the Cauchy residue 
theorem 11^. I n particular, it can be shown that under typical 
conditions llo . each b^^ (w, t) is analytic in the upper-half w- 
plane with b^^^j —0 as |aj| —>• oo, which allows us to evaluate 
Eq. ([TOb and obtain 

RuJoi^) ~ ^ ( 11 ) 


where b^^{t) = b^^{u},t)\^=u,^+iA^^- By setting uj = ujq+ 
in Eq. @, we finally obtain 



FIG. 2. (Color online) Bifurcation diagrams and Lyapunov expo¬ 
nent. Bifurcation diagrams of the (a) full and (b) reduced dynamics 
calculated by plotting the values of R{t) [for (a)] and B{t) [for (b)] 
evaluated at the times of surfaces of section piercing, (c) The largest 
Lyapunov exponent Alle as a function of K calculated using the re¬ 
duced system. 


£lrE<^o - uJo)/N]du} is nearly one. Replacing 

the quantity in ([T2l i hy bi {I = 1,..., iV) and regarding 
bi as representing the collective dynamics associated with os¬ 
cillators whose target frequencies fall in bin /, we achieve our 
dimensionality reduction. As we will see, N can be made 
much smaller than N, thus greatly reducing the computational 
complexity. To evaluate the degree of synchronization in the 
reduced system, we use the order parameter 


dbu. 


dt 


K 

- bi,R:, 

■ (12) 

«« - M 





UJ0,^0 


field version of the full system. Importantly, this formal¬ 
ism can be used to reduce the dimensionality of the system. 
Eor example. Ref. ll^ dealt with the effects of degree as¬ 
sortativity in the absence of frequency assortativity and used 
an equation analogous to (ITSl i to achieve dimensionality re¬ 
duction (i.e., Afc, bujo -t bk, Pu,o Pk, and 

o-k'^k)- Here we use Eq. (IT^ to investigate the 
effects of frequency assortativity in the network model de¬ 
scribed above, which has constant node degrees. We then 
use ([T2I 1 to achieve dimensionality reduction from the origi¬ 
nal N differential equations [Eq. ([T]i] to a much smaller num¬ 
ber N, by dividing the interval [wminjWmax] into N bins of 
width (wmax — Wmin)/-^, where the center frequency of the 
Z* bin is ujQ = co[, and Wmin and w^ax are chosen so that 


(13) 


which is the reduced-system analog to the order parameter de¬ 
fined in Eq. Q. Einally, we note that the distribution and 
assortativity function can be either constructed to rep¬ 

resent an ensemble of networks or sampled from a particular 
network realization, as we do here. 

Returning to the networks obtained by the model described 
above, we construct the corresponding reduced systems using 
N = 20 - a number small enough to significantly reduce the 
computational cost, but large enough to retain the dynamical 
complexity. Solving Eq. ([T2l i as K is increased from zero as 
in the full system, we plot B vs K for the non-assortative and 
assortative networks in Eigs. [TJa) and (b), respectively, using 
dashed red curves. We note that there is good agreement with 
the full system in both cases, and the reduced dynamics do 
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FIG. 3. (Color online) Chaos V5 periodicity. Time delay embeddings of the full system [R{t), R{t — r)] for (a) K = 0.09 and (b) 0.105 and 
the reduced system [B{t), B{t — r)] for (c) K — 0.09 and (d) 0.105 for r = 1. Also for K = 0.09 and 0.105, the dynamic correlations cfj 
[(e) and (f), respectively] as calculated from the full system. 


a particularly good job of reproducing the irregular oscilla¬ 
tions of the assortative network. Note that for small K the 
solid blue curve in Fig. [Ha) undergoes small fluctuations not 
present in the reduced mean-field solution (red dashed curve). 
These fluctuations become smaller (not shown) as N is in¬ 
creased keeping k/N fixed and can thus be explained as being 
due to finite network size iH. The irregular oscillations for 
K < 0.12 in Fig.lHb) turn out to be indicative of macroscopic 
chaos, as we will discuss below. 

Numerical investigations of chaos. We begin by construct¬ 
ing bifurcation diagrams of both the full and reduced system 
for the assortative case. To do so, we consider the time-delay 
embeddings {x,y) = [R{t),R{t — t)] and [B{t),B{t — r)]. 
For a given value of K, we record all the values of x when 
the line x = y is traversed from y > x to y < x after dis¬ 
carding transients. We use a value of t = 0.2, which is large 
enough to overcome small finite size fluctuations present in 
the full system, and small enough to capture the macroscopic 
features of the dynamics. We present the results in Fig. |2l 
plotting the bifurcation diagram of the full and reduced sys¬ 
tems in panels (a) and (b), respectively. Overall the results 
agree well, both indicating complex oscillations and intricate 
behavior leading up to transitions to periodic and then sta¬ 
tionary behavior. We note that the reduced model has thin 
regions of periodic behavior that we do not observe in the full 
system. We believe that this difference between Figs. |2|a) 
and |2jb) is due to the finite size induced noise-like fluctua¬ 
tions present in the real network but not in the reduced net¬ 
work (e.g., as also present for small K in Fig.[T]) and that this 
noise destroys the windows of periodicity seen in Fig. |2|b) 
(see ifl^ for a related finite network size noise phenomenon). 
In order to test this, we first add noise to the right-hand side 
of (fTTIi . which we then insert into (fTSI) . Simulations of this 
noisy model (not shown) confirm that even rather small noise 
is sufficient to destroy the thin regions of periodic behav¬ 
ior, while making a negligible effect on the dynamics for 
K > 0.12. Next, we take advantage of the lower complex¬ 
ity of the reduced system to calculate the largest Lyapunov 


exponent Xlle S, and plot the results in Fig. I2c). The 
largest Lyapunov exponent indicates that the system quickly 
transitions to chaos at a small coupling strength, and then in¬ 
termittently transitions between chaotic {Xlle > 0) and pe¬ 
riod {Xlle = 0) behavior. We also investigated the behavior 
of the Lyapunov dimension Dl by computing the whole Lya¬ 
punov spectrum (ordered Xlle = Ai > A 2 > ...), giving 
Dl = k + (Ai -F ■ • ■ -F Afe)/|Afc+iLwhere k is the largest 
index such that Ai -F ■ • • + Afc > 0 11^ . We And that Dl is 
large near the middle of the chaotic regime and significantly 
decreases as K is increased to approach the periodic regime 
(e.g., Dl ~ 13.64 and 3.71 at K = 0.05 and 0.095, respec¬ 
tively). 

Finally, we investigate the genesis of chaotic dynamics in 
assortative networks. To visualize and study the dynamics we 
consider time-delay embeddings and frequency correlations 
between pairs of oscillators, defined as c“ = (1 — |a;|® — 
uj‘f^\/\u!i — where we denote the effective frequency of 

oscillator i as = T~^ {t)dt for large enough 

to and T. In particular, c“ quantifies the degree to which os¬ 
cillators i and j evolve on their own {cf = 0) or in unison 
(c“ = 1). We choose examples of chaotic and periodic dy¬ 
namics that occur at K = 0.09 and 0.105, respectively, and 
plot in Fig. [3 the time-delay embeddings using r = 1 for 
the full system [left column, panels (a) and (b)] and for the 
reduced system [middle column, panels (c) and (d)]. Finally, 
we plot the frequency correlations calculated from the full sys¬ 
tems in the right column for both K = 0.9 (e) and K = 0.105 
(f), with cf = 0 and 1 corresponding to white and blue, 
respectively. The correlations are plotted so that the indices 

j increase with each oscillator’s target frequency. First, we 
note that the time-delay embeddings of the full [Figs. [3a) and 
[3b)] and reduced [Figs. [3[c) and[3[d)] dynamics match ex¬ 
tremely well for both the chaotic and periodic examples. Sec¬ 
ond, using the dynamic correlations we observe the formation 
of three [Figs. [3[e)] and two [Figs. [3[f)] large groups. We 
view such groups as meta-oscillators, and we interpret the ob¬ 
served dynamics as resulting from interactions of these meta 
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oscillators. When two (one) such meta-oscillators are present 
the macroscopic dynamics of the order parameter is observed 
to be periodic (steady), while chaos can (and typically does) 
occur when there are three or more groups. 

Discussion. In this Rapid Communication we have stud¬ 
ied the synchronization of assortative coupled oscillator net¬ 
works. Our main results are twofold. First, we have stud¬ 
ied frequency assortativity and found that this effect can in¬ 
duce large and robust regions of chaotic dynamics. We have 
supported our results using numerical simulations of regular 
graphs with constant degree in order to emphasize the impor¬ 
tance of frequency assortativity. Second, we showed that the 
dimensionality reduction method first presented in Ref. El 
can be extended to study this interesting case. We empha¬ 
size the strong correspondence between the dynamics of the 
full system and its low dimensional system reduction. In 
both contexts we have investigated the complicated dynam¬ 
ics that emerge using a combination of bifurcation diagrams, 
Lyapunov exponents, and time-delay embeddings. Finally, we 
discussed the genesis of chaos and showed that several locally 
synchronized groups; “meta-oscillators,” emerge in assorta¬ 
tive networks, allowing for chaos. 


Chaos in the macroscopic dynamics of networks of coupled 
phase oscillators has been observed previously, but in different 
contexts. These situation include one-way coupling between 
two groups of coupled oscillators dH, globally-coupled os¬ 
cillators with bimodal frequencies and an oscillating coupling 
strength oscillated in time and interacting communities 
of oscillators with different natural frequencies ll28l] . Our re¬ 
sults show that in very simple coupled oscillator networks 
with fixed parameters and no external driving, chaos can be 
induced merely by frequency assortativity. We attribute these 
chaotic dynamics to the formation of three or more meta¬ 
oscillators, which is in contrast to the periodic (often called 
standing-wave) behavior that emerges as the result of two 
meta-oscillators 


ACKNOWLEDGMENTS 

This work was supported by the James S. McDonnell Foun¬ 
dation (PSS) and ARO Grant W911NF1210101 (EO). 


[1] S. H. Strogatz, Sync: the Emerging Science of Spontaneous Or¬ 
der (Hypemion, 2003). 

[2] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, 
Phys. Rep. 469, 93 (2008). 

[3] J. Buck, Q. Rev. Biol. 63, 265 (1988). 

[4] L. Glass and M. C. Mackey, From Clocks to Chaos: The 
Rhythms of Life (Princeton University Press, Princeton, 1988). 

[5] A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, Nat. 
Phys. 9, 191 (2013). 

[6] S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and 
E. Ott, Nature (London) 438, 43 (2005). 

[7] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence 
(Springer, New York, 1984). 

[8] T. Ichinomiya, Phys. Rev. E 70, 026116 (2004). 

[9] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. E 71, 036151 
(2005). 

[10] E. Oh, K. Rho, H. Hong, and B. Kahng, Phys. Rev. E 72, 
047101 (2005). 

[11] A. Arenas, A. Diaz-Guilera, andC. J. Perez-Vicente, Phys. Rev. 
Lett. 96, 114102 (2006). 

[12] G. Bariev, T. M. Antonsen, and E. Ott, Chaos 21, 025103 

( 2011 ). 

[13] E. Ott and T. M. Antonsen, Chaos 18, 037113 (2008). 

[14] E. Ott and T. M. Antonsen, Chaos 19, 023117 (2009). 

[15] D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, 
Phys. Rev. Lett. 101, 084103 (2008). 

[16] L. M. Childs and S. H. Strogatz, Chaos 18, 043128 (2008). 

[17] E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So, and T. 
M. Antonsen, Phys. Rev. E 79, 026204 (2009). 


[18] W. S. Lee, E. Ott, and T. M. Antonsen, Phys. Rev. Lett. 103, 
044101 (2009) 

[19] P. S. Skardal, E. Ott, and J. G. Restrepo, Phys. Rev. E 84, 
036208 (2011). 

[20] E. Barreto, B. Hunt, E. Ott, and P. So, Phys. Rev. E 77, 036107 
(2008). 

[21] P. S. Skardal and J. G. Restrepo, Phys. Rev. E 85, 016208 

( 2012 ). 

[22] J. G. Restrepo and E. Ott, Europhys. Lett. 107, 60006 (2014). 

[23] M. E. J. Newman, Phys. Rev. Lett. 89, 208701 (2002). 

[24] M. E. J. Newman, Phys. Rev. E 67, 026126 (2003). 

[25] J. J. Pfeiffer III et ah. In Proceedings of the 23rd International 
World Wide Web Conference, 2014. 

[26] L. M. Alonso and G. B. Mindlin, Chaos 21, 023102 (2011). 

[27] P. So and E. Barreto, Chaos 21, 033127 (2011). 

[28] M. Komarov and A. Pikovsky, Phys. Rev. Lett. 110, 134101 
(2013). 

[29] M. Molloy and B. Reed, Random Struct. Algor. 6, 161 (1995). 

[30] E. Ereitag and R. Busam, Complex Analysis (Springer-Verlag, 
Berlin, 2009). 

[31] E. J. Hildebrand, M. A. Buice, and C. C. Chow, Phys. Rev. Lett. 
98, 054101 (2007). 

[32] E. Ott, Chaos in Dynamical Systems (Cambridge University 
Press, Cambridge, 2002). 

[33] J. L. Kaplan and J. A. Yorke, in Functional Differential Equa¬ 
tions and Approximations of Fixed Points edited by H.-O. Peit- 
gen and H.-O. Walther (Berlin: Springer-Verlag, Berlin, 1979) 
Volume 730, p. 204. 

[34] J. Gomez-Gardenes, S. Gomez, A. Arenas, and Y. Moreno, 
Phys. Rev. Lett. 106, 128701 (2011). 

[35] J. D. Crawford, J. Stat. Phys. 74, 1047 (1994). 



